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ABSTRACT 
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I Planets orbiting a planetesimal circumstellar disc can migrate inward from their ini- 
tial positions because of dynamical friction between planets and planetesimals. The 

'"^^1 migration rate depends on the disc mass and on its time evolution. Planets that are 

O ' embedded in long-lived planetesimal discs, having total mass of 10^'* — O.OlAfo, can 

^ I migrate inward a large distance and can survive only if the inner disc is truncated 



or because of tidal interaction with the star. In this case the semi-major axis, a, of 



5-H ' the planetary orbit is less than 0.1 AU. Orbits with larger a are obtained for smaller 

value of the disc mass or for a rapid evolution (depletion) of the disc. This model may 
explain several of the orbital features of the giant planets that were discovered in last 
years orbiting nearby stars as well as the metallicity enhancement found in several 
stars associated with short-period planets. 
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1 INTRODUCTION 

According to the most popular theory on the formation of giant planets in the solar system, planets were formed by accu- 
mulation of solid cores (Safronov 1969; Wetherill & Stewart 1989; Aarseth et al. 1993), known as planetesimals, in a gaseous 
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disc centered around the sun. When the core mass increases above IQMq, it begins a rapid accretion phase (Mizuno 1980; 
Bodenheimer & Pollack 1986) in which the protoplanet can capture a gas envelope Irom the protoplanetary disc leading to 
the formation of a giant planet (Pollack et al. 1996). Jupiter-mass planets may require most of the lifetime of the disc to 
accrete (lO'^-lO^ yr) (Zuckerman et al. 1995; Pollack et al. 1996). 

Protostellar discs around young stellar objects that have properties similar to that supposed for the solar nebula are com- 
mon: between 25 to 75% of young stellar objects in the Orion nebula seem to have discs (Prosser et al. 1994; McCaughrean 
& Stauffer 1994) with mass 10"^Mq < Afd < 10"^A'/q and size 40 ± 20 AU (Beckwith & Sargent 1996). Moreover recently 
several planetary companion orbiting extra-solar stars were discovered. The extrasolar planets census, updated at October 
2000, gives 58 planets. For reason of space, we report only a small number of them: the companions orbiting 51 Peg (Mayor 
& Queioz 1995), r Boo(Marcy et al. 1997 - San Francisco University Team hereafter SFSU), v And (SFSU), Cnc (SFSU), 
p CrB (Noyes et al. 1997 - AFOE team), HD 114762, 70 Vir, 16 Cyg, 47 UMa (Butler & Marcy 1996). In the above list, with 
the exception of 47 Uma, the new planets are all at distances < lAU. In recent years has been discovered other extrasolar 
planets orbiting at distances > lAU from the central star: 

e Eridani, HD210277, HD 82943, 14 Her, HD 190228, HD 222582, HD10697, HD 29587, representing only 15% of aU the 
planets. Three planets (51 Peg, r Boo, v And) are in extremely tight circular orbits with periods of a few days, two planets 
(p^ Cnc and p CrB)have circular orbits with periods of order tens of days and three planets with wider orbits (16 Cyg B, 70 
Vir and HD 114762) have very large eccentricities. The properties of these planets, most of which are Jupiter-mass objects, 
are difficult to explain using the quoted standard model for planet formation (Lissauer 1993; Boss 1995). This standard model 
predicts nearly circular planetary orbits, and giant planets distances > 1 AU from the central star, distance at which the 
temperature in the protostellar nebula is low enough for icy materials to condense (Boss 1995, 1996; Wuchterl 1993, 1996). 
Standard disc models show that at 0.05 AU, the temperature is about 2,000 K, which is too hot for the existence of any 
small solid particles. Moreover the ice condensation radius does not depend strongly on stellar mass, so that it does not move 
inward rapidly as the stellar mass decreases. For star masses, M« = 1, 0.5, 0.1 M©, the ice condensation radius moves inward 
from ~ 6 to ~ 4.5 AU(Boss 1995). Another problem with the in situ formation of a planetary companion is that even though 
the present evaporation rate is negligible, this effect would have been of major importance in the past. In fact during the 
early history of a planet, its radius was a factor ten larger than the present value, implying that the escape speed was much 
less than its present value. Hence evaporation mechanisms and ablation by the stellar wind might prevent its formation. The 
question that arises is: if such massive planets cannot form at the actual locations, how did they reach their actual position 
? Four mechanisms have been proposed to explain the quoted dilemma. 
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The first mechanism consists of a secular interaction with a distant binary companion (Holman et aL 1997; Mazeh et aL 
1996; Kiseleva & Eggleton 1997; Eggleton & Kiseleva 1997). While this mechanism can also produce significant eccentricities 
for the longer period extrasolar planets it is unable to explain objects like 51 Peg. In fact, 51 Peg has been extensively searched 
for a binary companion (Marcy et al. 1997), but none has been found. Consequently, in the particular case of 51 Peg, this 
mechanism is not responsible for the orbital decay. 

The second possible mechanism proposed to explain short period planets is dissipation in the protostellar nebula. Tidal 
interaction between a massive planet and a circumstellar disc gives rise to an angular momentum transfer between the disc 
and the planet (Goldreich & Tremaine 1979, 1980; Ward 1986; Lin et al. 1996; Ward 1997). The planet's motion in the disc 
excites density waves both interior and exterior to the planet. A torque originates from the attraction of the protoplanet for 
these non-axisymmetric density perturbations (Goldreich & Tremaine 1980). Density wave torques repel material on either 
side of the protoplanet's orbit and attempt to open a gap in the disc, whose size depends on the viscosity of the disc and 
inversely on the planet's mass (Lin & Papaloizou 1986; Takeuchi et al. 1986) (note that exists a minimum mass for gap 
opening, which is of the order of magnitude of Jovian planets mass, which prevents the nonsense of an infinite large gap for a 
zero- mass planet). If gap formation is successful (for example in the case of a Jupiter mass planet), the protoplanet becomes 
locked to the disc and must ultimately share its fate (Ward 1982; Lin & Papaloizou 1986, 1993). This mechanism is called 
type II drift. The situation is different if the object is not yet large enough to open and substain a gap. Also in this case the 
protoplanet migrates inwards but with a time-scale even smaller than that of type II drift (Ward 1997). This is called type I 
drift. In both cases, the rate of radial mobility of the planet, with respect to the central star, is indicated with the term 'drift 
velocity' (see Ward 1997) (in some cases, see the case of Neptun below, the drift velocity can be directed outwards). Since, 
in this model, the time-scale of migration is ~ lO'^lp-yr (Ward 1997), the migration has to switch off at a critical moment, 
if the planet has to stop close to the star without falling in it. The movement of the planet might be halted by short-range 
tidal or magnetic effects from the central star (Lin et al. 1996) (in any case, as shown by Murray et al. (1998), it is difficult 
to explain, by means of these stopping mechanisms, planets with semi- major axes a > 0.2 AU). 

A resonant interaction with a disc of planetesimals is another possible source of orbital migration. In this model, planet 
migration starts when the surface density of planetesimals, E satysfies the condition E > Ec, being Ec a critical value for Ec. 
The advantage of this mechanism is that the migration is halted naturally at short distances when the majority of perturbed 
planetesimals collide with the star. Moreover wide eccentric orbits can also be produced for planets more massive than ~ 3Mj. 
However the model has some disavantages, since the protoplanetary disc mass required for the migration of a Jupiter-mass 
planet to a ~ 0.1 AU is very large (Ford et al. 1999). 
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Finally, the fourth and last mechanism deals with dynamical instabilities in a system of giant planets (Rasio & Ford 1996). 
The orbits of planets could become unstable if the orbital radii evolve secularly at different rates or if the masses increase 
significantly as the planets accrete their gaseous envelopes (Lissauer 1993). In this model, the gravitational interaction between 
two planets, during evolution, (Gladman 1993; Chambers et al. 1996) can give rise to the ejection of one planet, leaving the 
other in eccentric orbit. The orbit of this can then circularize at an orbital separation of a few stellar radii (Rasio et al. 1996) 
if the inner planet has a sufhciently small pericenter distance. Simulations of multiple giant planets systems, showed that 
successive mergers between two or more planets can lead to the formation of a massive (> lOAfj) object in a wide eccentric 
orbit (Lin & Ida 1997). While it is almost certain that this mechanism operates in many systems with multiple planets, it is 
not clear if they can reproduce the fraction of systems similar to 51 Peg observed. 

In this paper we propose another model to explain the orbital parameters of extraterrestrial planets. The model is based on 
dynamical friction of a planet with a planetesimal disc. While the role of dynamical friction on the planetary accumulation 
process has been studied in several papers (Stewart & Kaula 1980; Horedt 1985; Stewart & Wetherill 1988), very few papers 
have studied its role on radial migration of planets or planetesimals in planetary discs (see Melita & Woolfson, 1996; Haghigh- 
ipour 1999). This attitude is due to the fact that, so far, many people have assumed a priori that radial migration due to 
dynamical friction is much slower than the damping of velocity dispersion due to dynamical friction. Therefore most studies on 
dynamical friction were concerned only with damping of velocity dispersion (damping of the eccentricity, e, and inclination, 
i), adopting local coordinates. Analytical works by Stewart & Wetherill (1988) and Ida (1990) adopted local coordinates. 
N-body simulation by Ida & Makino (1992) adopted non-local coordinates, but did not investigate radial migration. Moreover 
most models of the planetesimal disc assume that, except for the influence of aerodynamic drag, which loses its effectiveness 
for planetesimals larger than a few kilometers, the primary cause of radial migration is mutual scattering (Hayashi et al. 1977 
and Wetherill 1990). In order to calculate protoplanets migration, we apply the model introduced in Del Popolo et al. (1999) 
to study KBOs (Kuiper Belt Objects) migration, and we suppose that the gas in the disc is dissipated soon after the planet 
forms so that it has little effect on planet migration. We are particularly interested in studying the role of planetesimals in 
planet migration and the dependence of migration on the disc mass and on its evolution. 

The plan of the paper is the following: in Sect. 2 we introduce the model used to study radial migration. In Sect. 3 we 
show the assumption used in the simulation. In Sect. 4 we show the results that can be drawn from our calculations and 
finally the Sect. 5 is devoted to the conclusions. 
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2 PLANETS MIGRATION MODEL 

In a recent paper by Del Popolo et al. (1999), we studied how dynamical friction, due to small planetesimals, influences the 
evolution of KBOs having masses larger than 10^^ g. We found that the mean eccentricity of large mass particles is reduced 
by dynamical friction due to small mass particles in timescales shorter than the age of the solar system for objects of mass 
equal or larger than 10^'^ g. Moreover the dynamical drag, produced by dynamical friction of objects of masses > 10^* g, is 
responsible for the loss of angular momentum and the fall through more central regions in a timescale ~ lO' yr. Here we 
apply a similar model to study the planets radial migration. We suppose that a single planet moves in a planetesimal disc 
under the influence of the gravitational force of the Sun. The equation of motion of the planet can be written as: 



(Melita & Woolfson 1996), where the term F0 represents the force per unit mass from the Sun, while R is the dissipative 
force (the dynamical friction term-see Melita & Woolfson 1996). Calculations involving dynamical friction that are used to 
study planetesimal dynamics often use Chandrasekhar's theory (Stewart & Wheterill 1988; Ida 1990; Lissauer & Stewart 1992) 
for homogeneous and isotropic distribution of lighter particles. This choice is not the right one, since dynamical friction in 
discs differs from that in spherical isotropic three dimensional systems. This is because disc evolution is influenced by effects 
different than those producing the evolution of stellar systems: 

1) In a disc, the contribution to the friction coming from close encounters is comparable to that due to distant encounters 
(Donner & Sundelius 1993, Palmer et al. 1993). 

2) Collective effects in a disc are much stronger than those in a three-dimensional system (Thorne 1968). 

3) The peculiar velocities of planetesimals in a disc are small. This means that differential rotation of the disc dominates over 
planetesimals' relative velocities. 

4) The velocity dispersion of particles in a disc potential is anisotropic. 

We assume that the matter-distribution is disc-shaped and that it has a velocity distribution described by: 



(Hornung & al. 1985, Stewart & Wetherill 1988) where and (T|| are the velocity and the velocity dispersion in the direction 
parallel to the plane while v± and a± are those in the perpendicular direction. We suppose that (T|| and a± are constants and 
that their ratio is simply taken to be 2:1. Then according to Chandrasekhar (1968) and Binney (1977) we may write the force 



r = F© + R 



(1) 




(2) 
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components as: 



where 



2v27rnG Ig Amim2 (mi + 7712) 



2v27rnG Ig Amim2 (mi + m2 



(3) 
(4) 



B» 



exp 



exp 



1 



2cr2 1 + g 2o-2 1 - e2 + g 



1 



2cr2 1 + g 2cr2 1 - e2 + g 



[(l + g)2(l_e2+g)l/2] 

[(l + g)(l-e2 + Q)^/2] 



(5) 



(6) 



and 



e= (1- crx/o-||) 



(7) 



while n is the average spatial density, mi is the mass of the test particle, m2 is the mass of a field one, and log A is the 
Coulomb logarithm. The frictional drag on the test particles may be written as: 



F = — fc||Ui||e|| — k±Vi±ej_ 



(8) 



where ey and are two versors parallel and perpendicular to the disc plane. 

When B± > _B|| , the drag caused by dynamical friction will tend to increase the anisotropy of the velocity distribution 
of the test particles. In other words, the dynamical drag experienced by an object of mass mi moving through a less massive 
non-spherical distribution of objects of mass m2 is not directed in the direction of the relative motion (as in the case of 
spherically symmetric distribution of matter). Hence the already flat distribution of more massive objects will be further 
flattened during the evolution of the system (Binney 1977). As shown by Ida (1990), Ida & Makino (1992) and Del Popolo 
et al. (1999) damping of eccentricity and inclination is more rapid than radial migration so in this paper we deal only with 
radial migration and we assume that the planet has negligible inclination and eccentricity, ip ~ ep ~ and that the initial 
heliocentric distance of the planet is 5.2AU. The objects lying in the plane have no way of knowing that they are moving into 
a non-spherically symmetric potential. Hence we expect that the dynamical drag is directed in the direction opposite to the 
motion of the particle: 



F : 



-k\\v\\e\\ 



(9) 
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3 SIMULATION PARAMETERS 

In order to calculate the effect of dynamical friction on the orbital evolution of the planet, we suppose that a\^=2a± and that 
the dispersion velocities are constant. If the planetesimals attain dynamical equilibrium, their equilibrium velocity dispersion, 
Uni, would be comparable to the surface escape velocity of the dominant bodies (Safronov 1969) such that 

(^r (10) 



Jzn ' ^ i-'esc ' 



V Or, 

where 9 is the Safronov number, m« and r* are the mass and radius of the largest planetesimals, (note that the planetesimals 
velocity dispersion, am, now introduced, is the velocity dispersion to be used for calculating the a which is present in the 
dynamical friction force). If instead we consider a two-component system, consisting of one protoplanet and many equal- 
mass planetesimals the velocity dispersion of planetesimals in the neighborhood of the protoplanet depends on the mass of 
the protoplanet. When the mass of the planet, M, is < 10^^ g, the value of < >^^^ (being em the eccentricity of the 
planetesimals) is independent of Af therefore: 

em ~ 2O(2m/3M0)^/^ (11) 

(Ida & Makino 1993) where m is the mass of the planetesimals. When the mass of the planet reaches values larger than 
^q25_^q26 g ^ < e^ >^/^ is proportional to M^/^ 

em ~ 6(M/3M0)i/^ (12) 

(Ida & Makino 1993). As a consequence also the dispersion velocity in the disc is characterized by two regimes being it 
connected to the eccentricity by the equation: 

(Tm ~ (em + ii)^''^-Uc (13) 

where im is the inclination of planetesimals and Vc is the Keplerian circular velocity. The width of the heated region is roughly 
given by 4[(4/3)(em +im)a^ + l'2hlia^]^^^ (Ida & Makino 1993) where a is the semi-major axis and hm = i^M^y^'^ is the Hill 
radius of the protoplanet. The increase in velocity dispersion of planetesimals around the protoplanet decreases the dynamical 
friction force (see Eq. ^ and consequently increases the migration time-scale. 

In the simulation we assume that the planetesimals have all equal masses, m, and that m « M , M being the planet 
mass. This assumption does not affect the results, since dynamical friction does not depend on the individual masses of these 
particles but on their overall density. We also assume that the surface density in planetesimals varies as S(r) = Eq (lAU/r)''''^, 
where Eq, the surface density at 1 AU, is a free parameter. The total mass in the planetesimal disc within radius r is then: 
© 2000 RAS, MNRAS 000, §-?? 
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Md ~ 1.4 X lO-^(E0/lO=^g/cm2)(r-/AU)°-^ (14) 

We assume that ~ 1% of the disc mass is in the form of sohd particles (Stepinski & Valageas 1996). To be more precise, 
assuming metal abundance Z = 0.02, the disc mass in gas interior to Jupiter's orbit is 0.16Mq (So/lO^g/cm'^) (Murray et 
al. 1998). In our model the gas is almost totally dissipated when the planet begins to migrate and we assume that the value 
of disc mass reported in the following part of the paper is contained within 40 AU (Weidenshilling 1977). We integrated the 
equations of motion in heliocentric coordinates using the Bulirsch-Stoer method. 



4 RESULTS 

4.1 Migration in a non-evolving disc 

Our model starts with a fully formed gaseous giant planet of IMj at 5.2 AU. As mentioned in the previous section, the 
circumstellar disc is assumed to have a power-law radial density and to be axisymmetric. According to several evidences 
showing that the disc lifetimes range from 10^ yr to 10*^ yr (Strom et al. 1993; Ruden & Pollack 1991), we assume that the 
disc has a nominal effective lifetime of l(f years (Zuckerman et al. 1995). This assumption refers to the gas disc. Usually, 
this decline of gas mass near stars is more rapid than the decline in the mass of orbiting particulate matter (Zuckerman et 
al. 1995). Moreover the disc is populated by residual planetesimals for a longer period. We are interested in studying the 
migration due to interaction with planetesimals and for this reason we suppose that the gas is almost dissipated when the 
planet starts its migration. Since Jupiter-mass planets may require most of the lifetime of the disc to accrete (10^ to 10^ years) 
(and meanwhile the disc is subject to evolutionary changes) and since the disk is also subject to evolution after this time 
interval (Pollack et al. 1996; Zuckerman et al. 1995), we incorporated this possibility in our model by running also models 
allowing some disc to dissipate during the planets migration. We integrated the model introduced in the previous section for 
several values of the disc surface density or equivalently several disc masses: A/d = 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq. 

The results of this first set of calculations (assuming that the disc does not evolve) is shown in Fig. 1. The curves show 
the evolution of a 1 Mj planet in a disc with planetesimals surface density E = E0(lAU/r)^''^. The simulation is started with 
the planet at 5.2 AU and ip ~ Bp ~ 0. In any case, similarly to what showed by Murray et al. (1998), planets with masses 
M > 3Mj during their migration can increase the value of ep. The curves correspond, from bottom to top (short-dashed 
long-dashed line, dotted line, long dashed line, dot-short dashed line, solid line) to the following values of A/d: 0.01, 0.005, 
0.001, 0.0005, 0.0001 A/0. As expected the most massive disc (O.OlAf0) produces a rapid radial migration of the planet. 
Discs having masses lower than O.OIA/0 produce a smaller radial migration of the planet. In particular, we found that for 
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t(yr) 

Figure 1. The evolution of tlie a(t) of a Jupiter-mass planet, M = IMj in a disc with E = E0(lUA/r)'''''^ for several values of the 
planetesimal disc mass, Mj^ = O.OIMq (short-dashed long-dashed line), 0.005 Mq (dotted line), 0.001 Mq (long dashed line), 0.0005 
Mq (dot-short dashed line) and 0.0001 Mq (solid line). S is supposed to remain constant in time. 

Mu = O.OIMq the planet moves to 0.05 AU in ~ 4 X 10^ yr. The migration halts for the reason explained in the following 
(with the term 'halt' we mean that the planet has not had time to migrate any further, even if it is still migrating). If 
Md = 0.005,0.001,0.0005, O.OOOlAf© we have respectively, for the time needed to reach 0.05 AU: ~ 8 x 10^ yr, ~ 4 x 10** 
yr, ~ 7.5x — 10* yr, ~ 3.5 x 10^ yr. In other words, disc mass is one of the parameters that controls radial migration. An 
interesting feature of the model is that migration naturally halts without needing any peculiar mechanisms that avoid the 
planet from plunging into the central star. In fact as shown in Fig. 1 the migration time to reach ~ 0.05 AU increases with 
decreasing disk mass. If the disk mass is < O.OOOOSM© the time needed to reach the quoted position is larger than the age 
of the stellar system and the planet does not fall into the star. Then, the planet can halt its migration without falling in the 
star if the initial disc mass is < O.OOOOSM©. Even if the disc density does not fall below the critical value, the planet must 
halt at several J?, from the star surface (J?, is the stellar radius). In fact solid bodies cannot condense at distances < 7R*, 
and planetesimals cannot survive for a long time at distances < 2R* . When the planet arrives at this distance the dynamical 
friction force switches off and its migration stops. This means that the minimum value of the semi-major axis that a planet 
© 2000 RAS, MNRAS 000, §-?? 
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t(yr) 



Figure 2. The assumed evolution of tlie planetesimals mass. Tlie mass, from the initial value Mo, reduces to 10 ^Mo, in ~ 10^ yr (solid 
lino), ~ 10* yr (dotted line), ~ 3 X 10* (dashed line) 

can reach is ~ 0.03 AU. 



4.2 Migration in an evolving disc 



In the previous calculations, we supposed that the disc mass did not undergo time evolution but in reality the disc evolves as 
the planet moves inward and tends to dissipate. 

In our model, we are fundamentally interested in the evolution of solid matter and planetesimals in discs. To this aim, it is 
very important to note that the distribution of solid particles follows a global time evolution, which accompanies the time 
evolution of the gaseous component of the disc. Due to viscous torques, the gaseous disc spreads and its mass diminishes. If 
initially the solid particles are small and coupled to the gas, they decouple from it when they gain mass because of coagulation. 
(Stepinski & Valageas 1996). Particles having radius r < 0.1 cm, can be considered perfectly coupled to the gas, while those 
having r > 10^ cm can be considered completely decoupled from it and their mean velocities remain practically unchanged 
with time (Stepinski & Valageas 1996). Moreover, as shown in a recent study by Ida et al. (2000) the radial migration of a 
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Figure 3. Same as Fig.l but now we suppose that the disc mass decreases exponentially in 3 X 10* yr (see Fig. 2). 

planet of Jupiter-mass produce a very rapid capture of planetesimals in the 2:1 and 3:2 resonances: the resonance capture 
occurs if the migration time, Tmig, of the planet is Tmig > 10^ yr for 2:1 resonance and if Tmig > 10"^ yr for 3:2 resonance. If 
the result is correct this means that the disc should be rapidly depleted with a consequent rapid stopping of migration. As 
can be understood by what previously told, disc evolution depends on disc and system characteristics. 

In this paper we tried to take account of disc evolution by supposing that the total mass in solids decays with time from 
its original value to the present value as shown in Fig. 2 (see Stepinski & Valageas 1996). The results of this calculation are 
shown in Fig. 2 to Fig. 5. 

Fig. 2 shows the evolution of the disc mass used in the calculations of radial migration (see Stepinski & Valageas 1996, 
Fig 6). We suppose that the mass in the disc decreases exponentially with time from its original value A/o to 10~^Mo in 
3 X 10* yr (dashed line), 10* yr (dotted line) and 10^ yr (solid line). Fig. 3 is the same as Fig. 1 but now we suppose that the 
mass in the disc decreases as described. 

If Md = O.OIM©, the planet stops its migration at a ~ 0.1 AU, while if Mo = 0.005, 0.001, 0.0005, O.OOOIM© , we have 
a ~ 0.7,3.5,4.3,5 AU. Fig. 4 is obtained by supposing that the mass decreases, as previously quoted, in 10* yr. As can be 
© 2000 RAS, MNRAS 000, §-?? 
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Figure 4. Same as Fig. 2 but with planetesimals mass evolving in 10® yr (see Fig. 2) 



shown the planet embedded in a disc having M-q = 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq respectively migrates to 1.47, 2.7, 
4.58, 4.88, 5.1 AU. Finally in Fig. 5 the time-scale for disc evolution is 10^ yr. The planet embedded in a disc with Md ~ 0.01, 
0.005, 0.001, 0.0005, 0.0001 M© migrates respectively to 4.58, 4.88, 5.13, 5.17, 5.2 AU. 

Our distribution of final masses and heliocentric distances predicts that massive planets can be present at any heliocentric 
distances between their formation locations and extremely small orbits, depending on the initial mass of the disc and its 
evolution. As we shall show in the following, the model can explain the locations of not only the close companions (at < 0.1 
AU), but it can also reproduce other observed planets, including Jupiter. 



4.3 Comparison of the model results with observations 

Configuration of planets like r Bootis b, 51 Peg b, having very small semimajor axis, can be reproduced by models with no 
evolution (see Fig. 1) or high values of the initial disc mass and low time evolution (e.g., A/d ~ O.OIM©, tcvoi = 4 x 10* 
yr)(Boss 1996) and disc mass a bit lower can explain the configuration of planets like 55 Cnc b (a = 0.11 AU), p CrB b 
(a = 0.23 AU). The parameters of 47 UMa b (a = 2.11 AU) can be explained, for example, by supposing a low mass disc 
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Figure 5. Same as Fig. 4 but with planetesimals mass evolving in 10^ yr (see Fig. 2) 

and evolution (e.g., Md = 5 x 10~^Mq, fcvoi = 3 x 10® yr). Planets like 70 Vir b (a = 0.43, ep = 0.4; Msinip ~ 6.6Mj) 
and HD 114762 h {a — 0.3; ep = 0.25; Msinzp ~ lOMj) have high eccentricities and masses larger than Jupiter's. The low 
semimajor axis can be explained by radial migration, as shown, while the high value of eccentricities can be explained in a 
model of interaction planet- planetesimals like ours in a way similar to that shown by Murray et al. (1998). If a planet having 
mass A/ > 3Mj, which is the case of 70 Vir b and HD 114762, moves in a planetesimal disc during interactions, planetesimals 
scattered from their Hill sphere can be ejected with < 1, (where AE and AL are respectively the energy and angular 

momentum removed from a planet by the ejection of a planetesimal), and the eccentricity ep tends to increase. For sake of 
completeness we must say that there are also some systems a bit puzzling. For example v And has Af sinip = 0.68 and an high 
eccentricity, ep ~ 0.15, HD210277 has Afsinip — 1.36 and ep — 0.45. Such high eccentricities could be explained supposing 
an encounter with an object having mass ~ A/q or alternatively supposing that the systems are seen at small ip. In the case 
of the companion to 16 Cyg B (A/sinip — 1.5; ep ~ 0.68) the high eccentricity may be due to interactions with the stellar 
companion (Murray et al. 1998). 

Our Solar System could have been subject to giant planet migration. For example, a shrinkage of the orbit of Jupiter of 
© 2000 RAS, MNRAS 000, §-?? 
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Figure 6. Drift velocity, -jj, as a function of mass. Velocities are normalized to V = 2^^^^{^frQ. where Mq is an Earth mass, 
E the surface density, S7 is the angular velocity and a the dispersion velocity. The assumed conditions are those considered appropriate 
for the Jovian region and assuming that Md = O.OIMq. The line oc M correspond to model described by Ward (1997) and its behavior 
is valid till ~ O.lMg but beyond this value there is a transition to a behavior oc M" 

0.1-0.2 AU could naturally explain the depletion of the outer asteroid belt (Fernandez & Ip 1984; Lion & Malhotra 1997). 
Some of our model runs produce a 1 Mj planet that move from 5.2 AU inwards for a fraction of AU. The planetesimal disc 
enabling this small migration has a lifetime ~ 10^ yr, so that the disc gas must have disappeared soon after Jupiter fully 
formed (Boss 1996). 

This last result is obviously strictly valid only for a single planet orbiting around the Sun because in presence of several 
planets, migration becomes more complex. A close example is that of the solar system. In this case two planets, Uranus and 
Neptune were subject to outward migration, which is the opposite of what expected. Several models have been proposed to ex- 
plain this outward migration. A first model is connected to gravitational scattering between planet and residuals planetesimals 
(Malhotra 1993; Ida et al. 2000). A second model allowing Neptune outward migration is connected to the dissipation in the 
protostellar nebula. In this case both inwards and outwards planet migration are allowed. In fact in a viscous disc, gas inside a 
particular radius, known as the radius of maximum viscous stress, rmvs, drifts inwards as it loses angular momentum while gas 
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outside Tmvs expands outwards as it receives angular momentum (Lynden-Bell & Pringle 1974). Neptune's outwards migration 
is due to the fact that the gas in the Neptune forming region has a tendency to migrate outwards (Ruden & Lin 1986). In Fig. 
6 we show the drift velocity, as function of mass, Al. As shown, objects having masses < M® have velocity drift increasing 
as M, while after a threshold mass any further mass increases begins to slow down the drift. As the threshold is exceeded the 
motion fairly abruptly converts to a slower mode in which the drift velocity is independent of mass. As previously explained, 
this behaviour is due to the transition from a stage in which the dispersion velocity is independent of M to a stage in which 
it increases with M^^^ (Ida & Makino 1993). This last stage is known as the protoplanet-dominated stage. The phenomenon 
is equivalent to that predicted in the density wave approach (Goldreich & Tremaine 1980; Ward 1997). In this approach, the 
density wave torques repel material on either side of the protoplanet's orbit and attempt to open a gap in the disc. Only very 
large objects are able to open and sustain the gap. After gap formation, the drift rate of the planet is set by disc viscosity 
and is generally smaller than in absence of the gap. We stress that the decaying portion of the curve corresponding to the 
transition from the first to the second stage does not correspond to any particular model because following Ida & Makino 
(1993) we do not have information on the evolution of a in the transition regime. 

4.4 Enhancements of metallicity 

Another important point is that the dynamical processes leading to planets migration can also affect the evolution of the 
central star. Gonzales (1997, 1998a,b) showed that several stars with short-period planets have high metallicities, [Fe/H] > 0.2. 
Gonzales (1998b) proposes that their metallicities have been enhanced by the accretion of high Z-material which leads to the 
speculation that there may be a relationship between stars with higher metallicities and stars with planets. Alternatively, the 
correlation could arise, if metal-rich stars have metal-rich discs which are more likely to form planets. 

Several mechanisms have been proposed to explain the high metallicity of stars having extra-solar planets. One of this 
mechanism is related to the Lin et al. (1996) migration model but this model has two severe drawbacks. Models in which 
gas disc material accretes on the star are not able to significantly alter the observed metallicity because the disc has a 
metallicity slightly larger than the star. Good results are obtained in models in which asteroids or planetesimals accrete on to 
the star (Murray et al. 1998). As previously quoted, Murray et al. (1998) suggest that a giant planet can induce eccentricity 
growth among residual planetesimals through resonant interactions. Subsequent close encounters cause most of the affected 
planetesimals to be ejected outwards while the planet migrates inward. A substantial population of planetesimals could induce 
a Jupiter-mass planet to migrate a large distance inward. Neglecting any planetesimals that are scattered into the star until 
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the planet reaches its final orbit and assuming that the fraction of planetesimals scattered onto the star consist of only those 
that become planet crossing before colliding with the star, the mass accreted on the star is given by: 

Mace ~ .f{a)M/a (15) 
where 

/(a) ~ 0.15(a/AU)-°'^^'' (16) 

(Ford et al. 1999), where /(a) is the fraction of planetesimals scattered onto the star by the planet at distance a and a is 
a parameter ~ 0.5 — 1. As shown by Eq. |l^ the fraction of planetesimals scattered onto the star increases with decreasing 
distance from it. At a = 0.05 /(a) ^ 0.46 — 0.92 if a = 0.5 — 1. For example, the mass of planetesimals which would be 
scattered into 51 Peg, whose final orbit has a = 0.05 AU, is ISOAf® (Ford et al. 1999). By starting from a star with solar 
metallicity adding ISOM^ of asteroids to 51 Peg the observed [Fe/H] increases to 0.48 (Ford et al. 1999). The value found is 
an inferior limit because it takes account only of planetesimals scattered from the planet located in its final orbit. Moreover it 
is calculated for a value of Eq ~ 4000g/cm^ less than the maximum disc mass used in Murray et al. (1998). Using the largest 
values of Eq used by Murray et al. (1998) one expects a value of Afacc double than that previously quoted. 

Even if we assume Afacc ~ ISOA^q, this value is larger than that observed in 51 Peg ([Fe/H] ~ 0.21). Then following the 
Ford et al. (1999) model for 51 Peg, the prediction for metallicity abundance is larger than that observed. This means that 
unless most of the acquired heavy elements are able to difi'use in the radiative interior, the planetesimal-scattering scenario 
for orbital migration would require more than 90 % of the close encounters to result in the outward ejection of planetesimals 
(Sandquist et al. 1998). 

Our model predicts a smaller value for metallicity close to the observed value for 51 Peg. In fact, our model, similarly 
to that by Murray et al. (1998), explains the planets migration by planet-planetesimals interaction, but differently from the 
Murray's et al. (1998) our model needs less planetesimal mass for radial migration. When the planet reaches its final location, 
Eq. ^ ^0 (giving the quantity of planetesimals scattered in the star) can be applied after scaling Murray et al. (1998) result 
to reflect our disc mass. In the case of the more massive disc (A/d = O.OIA^q) we find Afacc ~ 40AfQ and [Fe/H] ~ 0.2. 

An important point to stress is that the plausibility of such an explanation depends, among other, on the size of the 
stellar convective envelope at the time of accretion. In order to be efficient, the accretion must take place sufficiently late in 
the stellar evolution when the outer convective envelope is shallow. Accretion taking place while the star was still on the pre- 
main-sequence, and consequently having a large convective envelope, would have little effect on stars' observed metallicities. 
The time it takes for a planet to migrate to a < 0.1 AU orbit, in our model, is > 10^ yr. Once the planet stops its migration. 
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planetesimals inside its orbit are quickly cleared out. As shown by Ford et al. (1999), the right time to produce the observed 
metallicities is at 2 — 3 x 10^ yr. This value can therefore be regarded as an important constraint for disc models. 

5 CONCLUSIONS 

The discovery (Mayor & Queloz 1995; Marcy & Butler 1996; Butler & Marcy 1996; Butler et al. 1997; Cochran et al. 1997; 
Noyes et al. 1997) of extra-solar planets has revitalized the discussion on the theory of planetary system formation and 
evolution. Although close giant planets formation may be theoretically possible (Wuchterl 1993; Wuchterl 1996), it requires 
the initial formation of a solid core of at least 5 -i- 10 M® which may be difficult to achieve very close to the parent star. It is 
therefore more likely that Jupiter-mass extra-solar planets cannot form at small heliocentric distances (Boss 1995; Guillot et 
al. 1996). After those discoveries, the idea that planets can migrate radially (Goldreich & Tremaine 1980; Ward & Hourigan 
1989; Lin & Papaloizou 1993; Lin et al. 1996) for long distances has been taken more seriously than was made in the past 
years. 

In this paper, we showed that dynamical friction between the planet and a planetesimals disc is an important mechanism 
for planet migration. We showed that migration of IMj planet to small heliocentric distances (0.05 AU) is possible for a disc 
with a total mass of 10^* 1O~^M0 (we remember that, according to Stepinski & Valageas 1996, and Murray et. al 1998, 
only ~ 1% of the disc mass is in the form of solid particles) if the planetesimal disc does not dissipate during the planet 
migration or if the disc has A/d > O.OIMq and the planetesimals are dissipated in ~ 10* yr. The model predicts that massive 
planets can be present at any heliocentric distances for the right value of disc mass and time evolution. 

We also showed that the drift velocity of planets and than the migration time are very similar to the predictions of the 
density wave approach (Ward 1997): the drift velocity increases as M for masses smaller than O.IM^b and is constant for larger 
masses. Finally we showed that the metallicity enhancement observed in several stars having extrasolar planets can also be 
explained, similarly to what proposed by Murray et al. (1998) and Ford et al. (1999) by means of scattering of planetesimals 
onto the parent star, after the planet reached its final configuration. Comparing our model with other models, that attempt 
to explain planets migration, we think our model has some advantages: 

1) differently from models based on the density wave theory (Goldreich & Tremaine 1980; Ward 1986, 1997), our model does 
not require a peculiar mechanism to stop the inward migration (Lin et al. 1996). Planet halt is naturally provided by the 
model. It can explain planets found at heliocentric distances of > 0.1 AU or planets having larger values of eccentricity. It 
can explain metallicity enhancements observed in stars having planets in short-period orbits; 
© 2000 RAS, MNRAS 000, §-?? 
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2) differently from Murray et al. (1998) model, our model shows that radial migration is possible with not too massive 
planetesimals disc (which is one of the drawbacks of the Murray's et al. (1998) model) and predicts the right metallicity 
enhancement. 
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ABSTRACT 

Planets orbiting a planetesimal circumstellar disc can migrate inward from their ini- 
tial positions because of dynamical friction between planets and planetesimals. The 
migration rate depends on the disc mass and on its time evolution. Planets that are 
embedded in long-lived planetesimal discs, having total mass of 10"'* — O.OIMq, can 
migrate inward a large distance and can survive only if the inner disc is truncated 
or because of tidal interaction with the star. In this case the semi-major axis, a, of 
the planetary orbit is less than O.IAU. Orbits with larger a are obtained for smaller 
value of the disc mass or for a rapid evolution (depletion) of the disc. This model may 
explain several of the orbital features of the giant planets that were discovered in last 
years orbiting nearby stars as well as the metallicity enhancement found in several 
stars associated with short-period planets. 

Key words: Planets and satellites: general; planetary system 



1 INTRODUCTION 



According to the most popular theory on the formation of giant planets in the solar system, planets were formed by accu- 



mulation of solid cores (Safronov 1969; Wetherill & Stewart 1989; Aarseth et al. 1993), known as planetesimals, in a gaseous 
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disc centered around the sun. When the core mass increases above IOMq, it begins a rapid accretion phase (Mizuno 1980; 
Bodenheimer & Pollack 1986) in which the protoplanet can capture a gas envelope from the protoplanetary disc leading to 
the formation of a giant planet (Pollack et al. 1996). Jupiter-mass planets may require most of the lifetime of the disc to 
accrete (lO'^-lO^ yr) (Zuckerman et al. 1995; Pollack et al. 1996). 

Protostellar discs around young stellar objects that have properties similar to that supposed for the solar nebula are com- 
mon: between 25 to 75% of young stellar objects in the Orion nebula seem to have discs (Prosser et al. 1994; McCaughrean 
& Stauffer 1994) with mass 10"^ < Md < Mq and size 40 ± 20 AU (Beckwith & Sargent 1996). Moreover recently 
several planetary companion orbiting extra-solar stars were discovered. The extrasolar planets census, updated at October 
2000, gives 58 planets. For reason of space, we report only a small number of them: the companions orbiting 51 Peg (Mayor 
& Queloz 1995), r Boo(Marcy et al. 1997 - San Francisco University Team hereafter SFSU), v And (SFSU), Cue (SFSU), 
p CrB (Noyes et al. 1997 - AFOE team), HD 114762, 70 Vir, 16 Cyg, 47 UMa (Butler & Marcy 1996). In the above list, with 
the exception of 47 Uma, the new planets are all at distances < lAU. In recent years has been discovered other extrasolar 
planets orbiting at distances > lAU from the central star: 

e Eridani, HD210277, HD 82943, 14 Her, HD 190228, HD 222582, HD10697, HD 29587, representing only 15% of all the 
planets. Three planets (51 Peg, r Boo, v And) are in extremely tight circular orbits with periods of a few days, two planets 
[p^ Cue and p CrB)have circular orbits with periods of order tens of days and three planets with wider orbits (16 Cyg B, 70 
Vir and HD 114762) have very large eccentricities. The properties of these planets, most of which are Jupiter-mass objects, 
are difficult to explain using the quoted standard model for planet formation (Lissauer 1993; Boss 1995). This standard model 
predicts nearly circular planetary orbits, and giant planets distances > 1 AU from the central star, distance at which the 
temperature in the protostellar nebula is low enough for icy materials to condense (Boss 1995, 1996; Wuchterl 1993, 1996). 
Standard disc models show that at 0.05 AU, the temperature is about 2,000 K, which is too hot for the existence of any 
small solid particles. Moreover the ice condensation radius does not depend strongly on stellar mass, so that it does not move 
inward rapidly as the stellar mass decreases. For star masses, M, = 1, 0.5, 0.1 Mq, the ice condensation radius moves inward 
from ~ 6 to ~ 4.5 AU(Boss 1995). Another problem with the in situ formation of a planetary companion is that even though 
the present evaporation rate is negligible, this effect would have been of major importance in the past. In fact during the 
early history of a planet, its radius was a factor ten larger than the present value, implying that the escape speed was much 
less than its present value. Hence evaporation mechanisms and ablation by the stellar wind might prevent its formation. The 
question that arises is: if such massive planets cannot form at the actual locations, how did they reach their actual position 
? Four mechanisms have been proposed to explain the quoted dilemma. 

© 2000 RAS, MNRAS 000, 1-?? 



Migration of giant planets in planetesimal discs. 3 

The first mechanism consists of a secuiar interaction with a distant binary companion (Hoiman et ai. 1997; Mazeh et ai. 
1996; Kiseieva & Eggieton 1997; Eggieton & Kiseieva 1997). Whiie this mechanism can aiso produce significant eccentricities 
for the longer period extrasolar planets it is unable to explain objects like 51 Peg. In fact, 51 Peg has been extensively searched 
for a binary companion (Marcy et al. 1997), but none has been found. Consequently, in the particular case of 51 Peg, this 
mechanism is not responsible for the orbital decay. 

The second possible mechanism proposed to explain short period planets is dissipation in the protostellar nebula. Tidal 
interaction between a massive planet and a circumstellar disc gives rise to an angular momentum transfer between the disc 
and the planet (Goldreich & Tremaine 1979, 1980; Ward 1986; Lin et al. 1996; Ward 1997). The planet's motion in the disc 
excites density waves both interior and exterior to the planet. A torque originates from the attraction of the protoplanet for 
these non-axisymmetric density perturbations (Goldreich & Tremaine 1980). Density wave torques repel material on either 
side of the protoplanet's orbit and attempt to open a gap in the disc, whose size depends on the viscosity of the disc and 
inversely on the planet's mass (Lin & Papaloizou 1986; Takeuchi et al. 1986) (note that exists a minimum mass for gap 
opening, which is of the order of magnitude of Jovian planets mass, which prevents the nonsense of an infinite large gap for a 
zero-mass planet). If gap formation is successful (for example in the case of a Jupiter mass planet), the protoplanet becomes 
locked to the disc and must ultimately share its fate (Ward 1982; Lin & Papaloizou 1986, 1993). This mechanism is called 
type 11 drift. The situation is different if the object is not yet large enough to open and substain a gap. Also in this case the 
protoplanet migrates inwards but with a time-scale even smaller than that of type 11 drift (Ward 1997). This is called type 1 
drift. In both cases, the rate of radial mobility of the planet, with respect to the central star, is indicated with the term 'drift 
velocity' (see Ward 1997) (in some cases, see the case of Neptun below, the drift velocity can be directed outwards). Since, 
in this model, the time-scale of mig ration is ~ 10^ ^yr (Ward 1997), the migration has to switch off at a critical moment, 
if the planet has to stop close to the star without falling in it. The movement of the planet might be halted by short-range 
tidal or magnetic effects from the central star (Lin et al. 1996) (in any case, as shown by Murray et al. (1998), it is difficult 
to explain, by means of these stopping mechanisms, planets with semi-major axes a > 0.2 AU). 

A resonant interaction with a disc of planetesimals is another possible source of orbital migration. In this model, planet 
migration starts when the surface density of planetesimals, E satysfies the condition E > Ec, being Ec a critical value for Ec. 
The advantage of this mechanism is that the migration is halted naturally at short distances when the majority of perturbed 
planetesimals collide with the star. Moreover wide eccentric orbits can also be produced for planets more massive than ~ 3Mj. 
However the model has some disavantages, since the protoplanetary disc mass required for the migration of a Jupiter-mass 
planet to a ~ 0.1 AU is very large (Ford et al. 1999). 
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Finally, the fourth and last mechanism deals with dynamical instabilities in a system of giant planets (Rasio & Ford 1996). 
The orbits of planets could become unstable if the orbital radii evolve secularly at different rates or if the masses increase 
significantly as the planets accrete their gaseous envelopes (Lissauer 1993). In this model, the gravitational interaction between 
two planets, during evolution, (Gladman 1993; Chambers et al. 1996) can give rise to the ejection of one planet, leaving the 
other in eccentric orbit. The orbit of this can then circularize at an orbital separation of a few stellar radii (Rasio et al. 1996) 
if the inner planet has a sufficiently small pericenter distance. Simulations of multiple giant planets systems, showed that 
successive mergers between two or more planets can lead to the formation of a massive (> lOMj) object in a wide eccentric 
orbit (Lin & Ida 1997). While it is almost certain that this mechanism operates in many systems with multiple planets, it is 
not clear if they can reproduce the fraction of systems similar to 51 Peg observed. 

In this paper we propose another model to explain the orbital parameters of extraterrestrial planets. The model is based on 
dynamical friction of a planet with a planetesimal disc. While the role of dynamical friction on the planetary accumulation 
process has been studied in several papers (Stewart & Kaula 1980; Horedt 1985; Stewart & Wetherill 1988), very few papers 
have studied its role on radial migration of planets or planetesimals in planetary discs (see Melita & Woolfson, 1996; Haghigh- 
ipour 1999). This attitude is due to the fact that, so far, many people have assumed a priori that radial migration due to 
dynamical friction is much slower than the damping of velocity dispersion due to dynamical friction. Therefore most studies on 
dynamical friction were concerned only with damping of velocity dispersion (damping of the eccentricity, e, and inclination, 
i), adopting local coordinates. Analytical works by Stewart & Wetherill (1988) and Ida (1990) adopted local coordinates. 
N-body simulation by Ida & Makino (1992) adopted non-local coordinates, but did not investigate radial migration. Moreover 
most models of the planetesimal disc assume that, except for the influence of aerodynamic drag, which loses its effectiveness 
for planetesimals larger than a few kilometers, the primary cause of radial migration is mutual scattering (Hayashi et al. 1977 
and Wetherill 1990). In order to calculate protoplanets migration, we apply the model introduced in Del Popolo et al. (1999) 
to study KBOs (Kuiper Belt Objects) migration, and we suppose that the gas in the disc is dissipated soon after the planet 
forms so that it has little effect on planet migration. We are particularly interested in studying the role of planetesimals in 
planet migration and the dependence of migration on the disc mass and on its evolution. 

The plan of the paper is the following: in Sect. 2 we introduce the model used to study radial migration. In Sect. 3 we 
show the assumption used in the simulation. In Sect. 4 we show the results that can be drawn from our calculations and 
finally the Sect. 5 is devoted to the conclusions. 
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2 PLANETS MIGRATION MODEL 

In a recent paper by Del Popolo et al. (1999), we studied how dynamical friction, due to small planetesimals, influences the 
evolution of KBOs having masses larger than 10^^ g. We found that the mean eccentricity of large mass particles is reduced 
by dynamical friction due to small mass particles in timescales shorter than the age of the solar system for objects of mass 
equal or larger than 10^'^ g. Moreover the dynamical drag, produced by dynamical friction of objects of masses > 10^* g, is 
responsible for the loss of angular momentum and the fall through more central regions in a timescale x 10^ yr. Here we 
apply a similar model to study the planets radial migration. We suppose that a single planet moves in a planetesimal disc 
under the influence of the gravitational force of the Sun. The equation of motion of the planet can be written as: 

r = F0 + R (1) 

(Melita & Woolfson 1996), where the term Fq represents the force per unit mass from the Sun, while R is the dissipative 
force (the dynamical friction term-see Melita & Woolfson 1996). Calculations involving dynamical friction that are used to 
study planetesimal dynamics often use Chandrasekhar's theory (Stewart & Wheterill 1988; Ida 1990; Lissauer & Stewart 1992) 
for homogeneous and isotropic distribution of lighter particles. This choice is not the right one, since dynamical friction in 
discs differs from that in spherical isotropic three dimensional systems. This is because disc evolution is influenced by effects 
different than those producing the evolution of stellar systems: 

1) In a disc, the contribution to the friction coming from close encounters is comparable to that due to distant encounters 
(Donner & Sundelius 1993, Palmer et al. 1993). 

2) Collective effects in a disc are much stronger than those in a three-dimensional system (Thorne 1968). 

3) The peculiar velocities of planetesimals in a disc are small. This means that differential rotation of the disc dominates over 
planetesimals' relative velocities. 

4) The velocity dispersion of particles in a disc potential is anisotropic. 

We assume that the matter-distribution is disc-shaped and that it has a velocity distribution described by: 



1 

ra(v,x) = ra(x) ( — 1 exp 



20-2 + 20-2 



(2) 



(Hornung & al. 1985, Stewart & Wetherill 1988) where ?)|| and (7|| are the velocity and the velocity dispersion in the direction 
parallel to the plane while v± and a± are those in the perpendicular direction. We suppose that (7|| and a± are constants and 
that their ratio is simply taken to be 2:1. Then according to Chandrasekhar (1968) and Binney (1977) we may write the force 
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components as: 



2v27rraG Ig Amim2 (mi + 77*2) 



where 



5,1 



and 



2v27rraG Ig Amim2 (mi + m2) 



(3) 



(4) 



exp 



exp 



2af.l + q 2af.l-e^+q 



dq 



2af.l + q 2af.l-e^+q 



[il + qfil-e^ + qy/^] 



dq 

[(l + g)(l-e2 + gf/^] 



(5) 



(6) 



2 / 2\0.5 



e = (1 - aj_/au) 



(7) 



while n is the average spatial density, mi is the mass of the test particle, m2 is the mass of a field one, and log A is the 
Coulomb logarithm. The frictional drag on the test particles may be written as: 



-A;||?)i||e|| - k±vi±e± 



(8) 



where ey and ej_ are two versors parallel and perpendicular to the disc plane. 

When B± > the drag caused by dynamical friction will tend to increase the anisotropy of the velocity distribution 
of the test particles. In other words, the dynamical drag experienced by an object of mass mi moving through a less massive 
non-spherical distribution of objects of mass m2 is not directed in the direction of the relative motion (as in the case of 
spherically symmetric distribution of matter). Hence the already flat distribution of more massive objects will be further 
flattened during the evolution of the system (Binney 1977). As shown by Ida (1990), Ida & Makino (1992) and Del Popolo 
et al. (1999) damping of eccentricity and inclination is more rapid than radial migration so in this paper we deal only with 
radial migration and we assume that the planet has negligible inclination and eccentricity, ip ~ ep ~ and that the initial 
heliocentric distance of the planet is 5.2AU. The objects lying in the plane have no way of knowing that they are moving into 
a non-spherically symmetric potential. Hence we expect that the dynamical drag is directed in the direction opposite to the 
motion of the particle: 



(9) 
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3 SIMULATION PARAMETERS 

In order to calculate the effect of dynamical friction on the orbital evolution of the planet, we suppose that (7||=2(7j_ and that 
the dispersion velocities are constant. If the planetesimals attain dynamical equilibrium, their equilibrium velocity dispersion, 
(7m, would be comparable to the surface escape velocity of the dominant bodies (Safronov 1969) such that 



(10) 



0r, 

where 9 is the Safronov number, m, and r, are the mass and radius of the largest planetesimals, (note that the planetesimals 
velocity dispersion, (7m, now introduced, is the velocity dispersion to be used for calculating the a which is present in the 
dynamical friction force). If instead we consider a two-component system, consisting of one protoplanet and many equal- 
mass planetesimals the velocity dispersion of planetesimals in the neighborhood of the protoplanet depends on the mass of 
the protoplanet. When the mass of the planet, M, is < 10^^ g, the value of < el, >^/^ (b eing em the eccentricity of the 
planetesimals) is independent of M therefore: 

em ~ 2O(2m/3M0)^^^ (11) 

(Ida & Makino 1993) where m is the mass of the planetesimals. When the mass of the planet reaches values larger than 
1025_io26 g at 1 AU, < el, >^l^ is proportional to M^^^: 

em ~6(M/3M0)'/' (12) 

(Ida & Makino 1993). As a consequence also the dispersion velocity in the disc is characterized by two regimes being it 
connected to the eccentricity by the equation: 

(7m ~ (em + I'D^^^fc (13) 

where im is the inclination of planetesimals and Vc is the Keplerian circular velocity. The width of the heated region is roughly 
given by 4[(4/3)(em + «m)a^ + (Ida & Makino 1993) where a is the semi-major axis and hyi = {^f^f^ is the Hill 

radius of the protoplanet. The increase in velocity dispersion of planetesimals around the protoplanet decreases the dynamical 
friction force (see Eq. 8) and consequently increases the migration time-scale. 

In the simulation we assume that the planetesimals have all equal masses, m, and that m << M, M being the planet 
mass. This assumption does not affect the results, since dynamical friction does not depend on the individual masses of these 
particles but on their overall density. We also assume that the surface density in planetesimals varies as E(r) = EQ(lAU/r)'^^^, 
where Eq, the surface density at 1 AU, is a free parameter. The total mass in the planetesimal disc within radius r is then: 
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Md ~ 1.4 X lO"^(E0/lO^g/cm^)(r/AU)''^^ (14) 

We assume that ~ 1% of the disc mass is in the form of solid particles (Stepinski & Valageas 1996). To be more precise, 
assuming metal abundance Z = 0.02, the disc mass in gas interior to Jupiter's orbit is 0.16MQ(EQ/10'^g/cm^) (Murray et 
al. 1998). In our model the gas is almost totally dissipated when the planet begins to migrate and we assume that the value 
of disc mass reported in the following part of the paper is contained within 40 AU (Weidenshilling 1977). We integrated the 
equations of motion in heliocentric coordinates using the Bulirsch-Stoer method. 



4 RESULTS 

4.1 Migration in a non-evolving disc 

Our model starts with a fully formed gaseous giant planet of IMj at 5.2 AU. As mentioned in the previous section, the 
circumstellar disc is assumed to have a power-law radial density and to be axisymmetric. According to several evidences 
showing that the disc lifetimes range from 10^ yr to 10'^ yr (Strom et al. 1993; Ruden & Pollack 1991), we assume that the 
disc has a nominal effective lifetime of 10^ years (Zuckerman et al. 1995). This assumption refers to the gas disc. Usually, 
this decline of gas mass near stars is more rapid than the decline in the mass of orbiting particulate matter (Zuckerman et 
al. 1995). Moreover the disc is populated by residual planetesimals for a longer period. We are interested in studying the 
migration due to interaction with planetesimals and for this reason we suppose that the gas is almost dissipated when the 
planet starts its migration. Since Jupiter-mass planets may require most of the lifetime of the disc to accrete (10^ to 10'^ years) 
(and meanwhile the disc is subject to evolutionary changes) and since the disk is also subject to evolution after this time 
interval (Pollack et al. 1996; Zuckerman et al. 1995), we incorporated this possibility in our model by running also models 
allowing some disc to dissipate during the planets migration. We integrated the model introduced in the previous section for 
several values of the disc surface density or equivalently several disc masses: Md = 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq. 

The results of this first set of calculations (assuming that the disc does not evolve) is shown in Fig. 1. The curves show 
the evolution of a 1 M j planet in a disc with planetesimals surface density E = EQ(lAU/r)'^^^. The simulation is started with 
the planet at 5.2 AU and ip ~ ep ~ 0. In any case, similarly to what showed by Murray et al. (1998), planets with masses 
M > 3Mj during their migration can increase the value of ep. The curves correspond, from bottom to top (short-dashed 
long-dashed line, dotted line, long dashed line, dot-short dashed line, solid line) to the following values of Md: 0.01, 0.005, 
0.001, 0.0005, 0.0001 Mq. As expected the most massive disc (O.OIMq) produces a rapid radial migration of the planet. 
Discs having masses lower than O.OIMq produce a smaller radial migration of the planet. In particular, we found that for 
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t(yr) 



Figure 1. The evolution of the a(t) of a Jupiter-mass planet, M = IMj in a disc with S = (lUA/r)'^/^ for several values of the 
planetesimal disc mass, M-q = O.OIMq (short-dashed long-dashed line), 0.005 Mq (dotted line), 0.001 Mq (long dashed line), 0.0005 
Mq (dot-short dashed line) and 0.0001 Mq (solid line). S is supposed to remain constant in time. 

Mt) = O.OIMq the planet moves to 0.05 AU in ~ 4 x 10'^ yr. The migration halts for the reason explained in the following 
(with the term 'halt' we mean that the planet has not had time to migrate any further, even if it is still migrating). If 
M-D = 0.005, 0.001, 0.0005, O.OOOIMq we have respectively, for the time needed to reach 0.05 AU: ~ 8 x 10^ yr, ~ 4 x 10** 
yr, ~ 7.5 X = 10^ yr, ~ 3.5 x 10^ yr. In other words, disc mass is one of the parameters that controls radial migration. An 
interesting feature of the model is that migration naturally halts without needing any peculiar mechanisms that avoid the 
planet from plunging into the central star. In fact as shown in Fig. 1 the migration time to reach ~ 0.05 AU increases with 
decreasing disk mass. If the disk mass is < O.OOOOSMq the time needed to reach the quoted position is larger than the age 
of the stellar system and the planet does not fall into the star. Then, the planet can halt its migration without falling in the 
star if the initial disc mass is < O.OOOOSMq. Even if the disc density does not fall below the critical value, the planet must 
halt at several i?, from the star surface (i?, is the stellar radius). In fact solid bodies cannot condense at distances < 7i?,, 
and planetesimals cannot survive for a long time at distances < 2i?,. When the planet arrives at this distance the dynamical 
friction force switches off and its migration stops. This means that the minimum value of the semi-major axis that a planet 
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t(yr) 



Figure 2. The assumed evolution of the planetesimals mass. The mass, from the initial value -Mo, reduces to 10 ^ Mo , in ~ 10^ yr (solid 
line), ~ 10* yr (dotted line), ~ 3 X 10* (dashed line) 



can reach is ~ 0.03 AU. 



4.2 Migration in an evolving disc 



In the previous calculations, we supposed that the disc mass did not undergo time evolution but in reality the disc evolves as 
the planet moves inward and tends to dissipate. 

In our model, we are fundamentally interested in the evolution of solid matter and planetesimals in discs. To this aim, it is 
very important to note that the distribution of solid particles follows a global time evolution, which accompanies the time 
evolution of the gaseous component of the disc. Due to viscous torques, the gaseous disc spreads and its mass diminishes. If 
initially the solid particles are small and coupled to the gas, they decouple from it when they gain mass because of coagulation. 
(Stepinski & Valageas 1996). Particles having radius r < 0.1 cm, can be considered perfectly coupled to the gas, while those 
having r > 10^ cm can be considered completely decoupled from it and their mean velocities remain practically unchanged 
with time (Stepinski & Valageas 1996). Moreover, as shown in a recent study by Ida et al. (2000) the radial migration of a 

© 2000 RAS, MNRAS 000, 1-?? 



Migration of giant planets in planetesimal discs. 11 



0.1 





- 
- 

- 


" - ^ --^ ^ 

\ 

\ 

\ 


- 
- 

- 




\ 












\ 




- 


\ 


- 




\ 






\ 






\ 






\ 






\ 












1 1 1 





5x10^ 5x10" 5xl0' 5x10" 

t(yr) 

Figure 3. Same as Fig.l but now we suppose that the disc mass decreases exponentially in 3 X 10^ yr (see Fig. 2). 

planet of Jupiter-mass produce a very rapid capture of planetesimals in the 2:1 and 3:2 resonances: the resonance capture 
occurs if the migration time, Tmig, of the planet is Tmig > 10* yr for 2:1 resonance and if Tmig > 10^ yr for 3:2 resonance. If 
the result is correct this means that the disc should be rapidly depleted with a consequent rapid stopping of migration. As 
can be understood by what previously told, disc evolution depends on disc and system characteristics. 

In this paper we tried to take account of disc evolution by supposing that the total mass in solids decays with time from 
its original value to the present value as shown in Fig. 2 (see Stepinski & Valageas 1996). The results of this calculation are 
shown in Fig. 2 to Fig. 5. 

Fig. 2 shows the evolution of the disc mass used in the calculations of radial migration (see Stepinski & Valageas 1996, 
Fig 6). We suppose that the mass in the disc decreases exponentially with time from its original value Mo to 10~'^Mo in 
3 X 10* yr (dashed line), 10* yr (dotted line) and 10'^ yr (solid line). Fig. 3 is the same as Fig. 1 but now we suppose that the 
mass in the disc decreases as described. 

If M-D = O.OIMq, the planet stops its migration at a ~ 0.1 AU, while if Mb = 0.005, 0.001, 0.0005, O.OOOIM© , we have 
a ~ 0.7,3.5,4.3,5 AU. Fig. 4 is obtained by supposing that the mass decreases, as previously quoted, in 10* yr. As can be 
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Figure 4. Same as Fig. 2 but with planetesimals mass evolving in 10^ yr (see Fig. 2) 

shown the planet embedded in a disc having Md = 0.01, 0.005, 0.001, 0.0005, 0.0001 Mq respectively migrates to 1.47, 2.7, 
4.58, 4.88, 5.1 AU. Finally in Fig. 5 the time-scale for disc evolution is 10'^ yr. The planet embedded in a disc with Md = 0.01, 
0.005, 0.001, 0.0005, 0.0001 Mq migrates respectively to 4.58, 4.88, 5.13, 5.17, 5.2 AU. 

Our distribution of final masses and heliocentric distances predicts that massive planets can be present at any heliocentric 
distances between their formation locations and extremely small orbits, depending on the initial mass of the disc and its 
evolution. As we shall show in the following, the model can explain the locations of not only the close companions (at < 0.1 
AU), but it can also reproduce other observed planets, including Jupiter. 

4.3 Comparison of the model results with observations 

Configuration of planets like r Bootis b, 51 Peg b, having very small semimajor axis, can be reproduced by models with no 
evolution (see Fig. 1) or high values of the initial disc mass and low time evolution (e.g., Md = O.OIMq, tevoi = 4 x 10* 
yr)(Boss 1996) and disc mass a bit lower can explain the configuration of planets like 55 Cue b (a = 0.11 AU), p CrB b 
(a = 0.23 AU). The parameters of 47 UMa b (a = 2.11 AU) can be explained, for example, by supposing a low mass disc 
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Figure 5. Same as Fig. 4 but with planetesimals mass evolving in 10^ yr (see Fig. 2) 

and evolution (e.g., Md = 5 x IO'^Mq, tevoi = 3 x 10** yr). Planets like 70 Vir b (a = 0.43, Cp = 0.4; M sin Zp ~ 6.6Mj) 
and HD 114762 b (a = 0.3; Cp = 0.25; M sin ip ~ lOMj) have high eccentricities and masses larger than Jupiter's. The low 
semimajor axis can be explained by radial migration, as shown, while the high value of eccentricities can be explained in a 
model of interaction planet-planetesimals like ours in a way similar to that shown by Murray et al. (1998). If a planet having 
mass M > 3Mj, which is the case of 70 Vir b and HD 114762, moves in a planetesimal disc during interactions, planetesimals 
scattered from their Hill sphere can be ejected with |-^^| < 1, (where AE and AL are respectively the energy and angular 
momentum removed from a planet by the ejection of a planetesimal), and the eccentricity Cp tends to increase. For sake of 
completeness we must say that there are also some systems a bit puzzling. For example v And has M sin ip = 0.68 and an high 
eccentricity, Cp ~ 0.15, HD210277 has M sin ip = 1.36 and Cp = 0.45. Such high eccentricities could be explained supposing 
an encounter with an object having mass ~ or alternatively supposing that the systems are seen at small ip. In the case 
of the companion to 16 Cyg B (M sin ip = 1.5; Cp ~ 0.68) the high eccentricity may be due to interactions with the stellar 
companion (Murray et al. 1998). 

Our Solar System could have been subject to giant planet migration. For example, a shrinkage of the orbit of Jupiter of 
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Figure 6. Drift velocity, as a function of mass. Velocities are normalized to V = 2^^^(^ frQ where M® is an Earth mass, 




S the surface density, Q is the angular velocity and (T the dispersion velocity. The assumed conditions are those considered appropriate 
for the Jovian region and assuming that M-q = O.OIMq. The line oc M correspond to model described by Ward (1997) and its behavior 
is valid till ~ O.IMq^ but beyond this value there is a transition to a behavior oc M'^ 

0.1-0.2 AU could naturally explain the depletion of the outer asteroid belt (Fernandez & Ip 1984; Lion & Malhotra 1997). 
Some of our model runs produce a 1 Mj planet that move from 5.2 AU inwards for a fraction of AU. The planetesimal disc 
enabling this small migration has a lifetime ~ 10'^ yr, so that the disc gas must have disappeared soon after Jupiter fully 
formed (Boss 1996). 

This last result is obviously strictly valid only for a single planet orbiting around the Sun because in presence of several 
planets, migration becomes more complex. A close example is that of the solar system. In this case two planets, Uranus and 
Neptune were subject to outward migration, which is the opposite of what expected. Several models have been proposed to ex- 
plain this outward migration. A first model is connected to gravitational scattering between planet and residuals planetesimals 
(Malhotra 1993; Ida et al. 2000). A second model allowing Neptune outward migration is connected to the dissipation in the 
protostellar nebula. In this case both inwards and outwards planet migration are allowed. In fact in a viscous disc, gas inside a 
particular radius, known as the radius of maximum viscous stress, rmvs, drifts inwards as it loses angular momentum while gas 
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outside rmvs expands outwards as it receives angular momentum (Lynden-Bell & Pringle 1974). Neptune's outwards migration 
is due to the fact that the gas in the Neptune forming region has a tendency to migrate outwards (Ruden & Lin 1986). In Fig. 
6 we show the drift velocity, as function of mass, M. As shown, objects having masses < Mg^ have velocity drift increasing 
as M , while after a threshold mass any further mass increases begins to slow down the drift. As the threshold is exceeded the 
motion fairly abruptly converts to a slower mode in which the drift velocity is independent of mass. As previously explained, 
this behaviour is due to the transition from a stage in which the dispersion velocity is independent of M to a stage in which 
it increases with M^'^ (Ida & Makino 1993). This last stage is known as the protoplanet-dominated stage. The phenomenon 
is equivalent to that predicted in the density wave approach (Goldreich & Tremaine 1980; Ward 1997). In this approach, the 
density wave torques repel material on either side of the protoplanet's orbit and attempt to open a gap in the disc. Only very 
large objects are able to open and sustain the gap. After gap formation, the drift rate of the planet is set by disc viscosity 
and is generally smaller than in absence of the gap. We stress that the decaying portion of the curve corresponding to the 
transition from the first to the second stage does not correspond to any particular model because following Ida & Makino 
(1993) we do not have information on the evolution of a in the transition regime. 

4.4 Enhancements of metallicity 

Another important point is that the dynamical processes leading to planets migration can also affect the evolution of the 
central star. Gonzales (1997, 1998a, b) showed that several stars with short-period planets have high metallicities, [Fe/H] > 0.2. 
Gonzales (1998b) proposes that their metallicities have been enhanced by the accretion of high Z-material which leads to the 
speculation that there may be a relationship between stars with higher metallicities and stars with planets. Alternatively, the 
correlation could arise, if metal-rich stars have metal-rich discs which are more likely to form planets. 

Several mechanisms have been proposed to explain the high metallicity of stars having extra-solar planets. One of this 
mechanism is related to the Lin et al. (1996) migration model but this model has two severe drawbacks. Models in which 
gas disc material accretes on the star are not able to significantly alter the observed metallicity because the disc has a 
metallicity slightly larger than the star. Good results are obtained in models in which asteroids or planetesimals accrete on to 
the star (Murray et al. 1998). As previously quoted, Murray et al. (1998) suggest that a giant planet can induce eccentricity 
growth among residual planetesimals through resonant interactions. Subsequent close encounters cause most of the affected 
planetesimals to be ejected outwards while the planet migrates inward. A substantial population of planetesimals could induce 
a Jupiter-mass planet to migrate a large distance inward. Neglecting any planetesimals that are scattered into the star until 
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the planet reaches its final orbit and assuming that the fraction of planetesimals scattered onto the star consist of only those 
that become planet crossing before colliding with the star, the mass accreted on the star is given by: 

Mace ~ f{a}M/a (15) 
where 

/(a) ~ 0.15(a/AU)-'' "'' (16) 

(Ford et al. 1999), where f{a) is the fraction of planetesimals scattered onto the star by the planet at distance a and a is 
a parameter ~ 0.5 — 1. As shown by Eq. 16 the fraction of planetesimals scattered onto the star increases with decreasing 
distance from it. At a = 0.05 f{a) ~ 0.46 — 0.92 if a = 0.5 — 1. For example, the mass of planetesimals which would be 
scattered into 51 Peg, whose final orbit has a = 0.05 AU, is 13OM0 (Ford et al. 1999). By starting from a star with solar 
metallicity adding 13OM0 of asteroids to 51 Peg the observed [Fe/H] increases to 0.48 (Ford et al. 1999). The value found is 
an inferior limit because it takes account only of planetesimals scattered from the planet located in its final orbit. Moreover it 
is calculated for a value of Eq ~ 4000g/cm^ less than the maximum disc mass used in Murray et al. (1998). Using the largest 
values of Eq used by Murray et al. (1998) one expects a value of Mace double than that previously quoted. 

Even if we assume Mace ~ 130Mq, this value is larger than that observed in 51 Peg ([Fe/H] ~ 0.21). Then following the 
Ford et al. (1999) model for 51 Peg, the prediction for metallicity abundance is larger than that observed. This means that 
unless most of the acquired heavy elements are able to diffuse in the radiative interior, the planetesimal-scattering scenario 
for orbital migration would require more than 90 % of the close encounters to result in the outward ejection of planetesimals 
(Sandquist et al. 1998). 

Our model predicts a smaller value for metallicity close to the observed value for 51 Peg. In fact, our model, similarly 
to that by Murray et al. (1998), explains the planets migration by planet-planetesimals interaction, but differently from the 
Murray's et al. (1998) our model needs less planetesimal mass for radial migration. When the planet reaches its final location, 
Eq. 15 -^16 (giving the quantity of planetesimals scattered in the star) can be applied after scaling Murray et al. (1998) result 
to reflect our disc mass. In the case of the more massive disc (Md = O.OIMq) we find Mace ~ 4OM0 and [Fe/H] ~ 0.2. 

An important point to stress is that the plausibility of such an explanation depends, among other, on the size of the 
stellar convective envelope at the time of accretion. In order to be efficient, the accretion must take place sufficiently late in 
the stellar evolution when the outer convective envelope is shallow. Accretion taking place while the star was still on the pre- 
main-sequence, and consequently having a large convective envelope, would have little effect on stars' observed metallicities. 
The time it takes for a planet to migrate to a < 0.1 AU orbit, in our model, is > 10'^ yr. Once the planet stops its migration, 
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planetesimals inside its orbit are quickly cleared out. As shown by Ford et al. (1999), the right time to produce the observed 
metallicities is at 2 — 3 x 10'^ yr. This value can therefore be regarded as an important constraint for disc models. 

5 CONCLUSIONS 

The discovery (Mayor & Queloz 1995; Marcy & Butler 1996; Butler & Marcy 1996; Butler et al. 1997; Cochran et al. 1997; 
Noyes et al. 1997) of extra-solar planets has revitalized the discussion on the theory of planetary system formation and 
evolution. Although close giant planets formation may be theoretically possible (Wuchterl 1993; Wuchterl 1996), it requires 
the initial formation of a solid core of at least 5-^10 which may be difficult to achieve very close to the parent star. It is 
therefore more likely that Jupiter-mass extra-solar planets cannot form at small heliocentric distances (Boss 1995; Guillot et 
al. 1996). After those discoveries, the idea that planets can migrate radially (Goldreich & Tremaine 1980; Ward & Hourigan 
1989; Lin & Papaloizou 1993; Lin et al. 1996) for long distances has been taken more seriously than was made in the past 
years. 

In this paper, we showed that dynamical friction between the planet and a planetesimals disc is an important mechanism 
for planet migration. We showed that migration of IMj planet to small heliocentric distances (0.05 AU) is possible for a disc 
with a total mass of 10~* 10~^Mq (we remember that, according to Stepinski & Valageas 1996, and Murray et. al 1998, 
only ~ 1% of the disc mass is in the form of solid particles) if the planetesimal disc does not dissipate during the planet 
migration or if the disc has Mb > O.OIMq and the planetesimals are dissipated in ~ 10* yr. The model predicts that massive 
planets can be present at any heliocentric distances for the right value of disc mass and time evolution. 

We also showed that the drift velocity of planets and than the migration time are very similar to the predictions of the 
density wave approach (Ward 1997): the drift velocity increases as M for masses smaller than O.IM0 and is constant for larger 
masses. Finally we showed that the metallicity enhancement observed in several stars having extrasolar planets can also be 
explained, similarly to what proposed by Murray et al. (1998) and Ford et al. (1999) by means of scattering of planetesimals 
onto the parent star, after the planet reached its final configuration. Comparing our model with other models, that attempt 
to explain planets migration, we think our model has some advantages: 

1) differently from models based on the density wave theory (Goldreich & Tremaine 1980; Ward 1986, 1997), our model does 
not require a peculiar mechanism to stop the inward migration (Lin et al. 1996). Planet halt is naturally provided by the 
model. It can explain planets found at heliocentric distances of > 0.1 AU or planets having larger values of eccentricity. It 
can explain metallicity enhancements observed in stars having planets in short-period orbits; 
© 2000 RAS, MNRAS 000, 1-?? 
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2) differently from Murray et al. (1998) model, our model shows that radial migration is possible with not too massive 
planetesimals disc (which is one of the drawbacks of the Murray's et al. (1998) model) and predicts the right metallicity 
enhancement. 
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